###########################################################################
# IMPORTANT: Run via '00_master.R' or ensure WD is 'Codes' folder
###########################################################################

# --- 1. Path Management (Replication Safe) ---
# If this script is run standalone (not from Master), define defaults.
if (!exists("figures_dir")) figures_dir <- "Figures"
if (!exists("data_dir"))    data_dir    <- "Data"

# Create Figures folder if running standalone and it's missing
if (!dir.exists(figures_dir)) dir.create(figures_dir, showWarnings = FALSE)

# --- 2. Dependency Check ---
# This script relies on objects created in the make.tables.R script.
if (!exists("experiment.df")) {
  stop("CRITICAL ERROR: Required data objects are missing. You must run '01_make.tables.R' before running this script.")
}

###################################
############ Main Text ############ 
###################################


## Figure 3: means and distributions of outcome variables ##

#function to create summary and plot
make_summary_plot <- function(data, yvar, title) {
  
  summary_df <- data %>% 
    group_by(conditions) %>%
    do(tidy(lm_robust(reformulate("1", yvar), data = .))) %>%
    mutate(yval = estimate, significant = p.value <= 0.05)
  
  ggplot(data, aes(x = conditions, y = .data[[yvar]])) +
    geom_point(position = position_jitter(width = .1, height = .1, seed = 5208), 
               alpha = 0.1, stroke = 0) +
    geom_point(data = summary_df, aes(y = yval), colour = "#d62728") +
    geom_errorbar(data = summary_df,
                  aes(y = yval, ymin = conf.low, ymax = conf.high),
                  width = 0, colour = "#d62728") +
    theme_bw() +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text = element_text(size = 14, face = "bold", margin = margin(b = 5)),
          axis.title.x = element_blank(),
          axis.text.y = element_text(size = 12, color = "black"), 
          axis.text.x = element_text(size = 12, color = "black"),
          plot.title = element_text(hjust = 0.5)) +
    ylab("") + 
    ggtitle(title)
}


plot_Y1 <- make_summary_plot(experiment.df, "vote_han_change", "Δ (Vote for Han)")
plot_Y2 <- make_summary_plot(experiment.df, "relative_thermo_change", "Δ (Candidate Evaluations)")
plot_Y3 <- make_summary_plot(experiment.df, "prc_index_change", "Δ (Pro-PRC Index)")
combine_plot_raw <- plot_grid(plot_Y1,plot_Y2,plot_Y3, ncol=3) #combine into one figure
combine_plot_raw
ggsave(file.path(figures_dir, "Figure 3.pdf"), combine_plot_raw, width = 14, height = 6, dpi = 300, device = my_pdf_device)


## Figure 4: Treatment Effects
# Function to process models for each outcome
make_main_plot <- function(existing_model, itt_model, min_model, full_model, title) {
  bind_rows(
    tidy(existing_model) %>% mutate(model = "Existing"),
    tidy(itt_model) %>% filter(term != "(Intercept)") %>% mutate(model = "ITT"),
    tidy(min_model) %>% filter(term != "(Intercept)") %>% mutate(model = "Minimal"),
    tidy(full_model) %>% filter(term != "(Intercept)") %>% mutate(model = "Full")
  ) %>% mutate(title = title)
}

# Process all outcomes
main.plot <- bind_rows(
  make_main_plot(vote.han.existing, vote.han.itt, vote.han.min, vote.han.full, "Δ (Vote for Han)"),
  make_main_plot(favhan.existing, favhan.itt, favhan.min, favhan.full, "Δ (Candidate Evaluations)"),
  make_main_plot(china.existing, china.itt, china.min, china.full, "Δ (Pro-PRC Index)"))

# Set outcome order
main.plot$title <- factor(main.plot$title, 
                          levels = c("Δ (Vote for Han)", 
                                     "Δ (Candidate Evaluations)",
                                     "Δ (Pro-PRC Index)"))
main.plot$model <- factor(main.plot$model, 
                          levels = c("Existing", "Full", "Minimal", "ITT"))

# Plot
ggplot(main.plot, aes(x = model, y = estimate)) +
  geom_hline(yintercept = 0, color = "gray60", linetype = 'dashed') +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  linewidth = 1,
                  colour = "#d62728") + 
  scale_x_discrete(labels = c("Existing CT\nConsumers", 
                              "Treatment-on-the\nTreated (Full)", 
                              "Treatment-on-the\nTreated (Minimum)", 
                              "Intent-to-Treat")) +
  facet_grid(. ~ title, scales = "free") + coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        legend.key.width = unit(2, "lines"),
        strip.background = element_blank(),
        strip.text = element_text(size = 14, face = "bold", margin = margin(b = 5)),
        axis.text.y = element_text(size = 12, color = "black"), 
        axis.text.x = element_text(size = 12, color = "black"), 
        axis.title.x = element_text(size = 12),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank()) +
  xlab("") + ylab("")
ggsave(file.path(figures_dir, "Figure 4.pdf"), width = 12, height = 4, dpi = 300, device = my_pdf_device)


## Figure 5: Treatment Effects by Political Attentiveness
make_interaction_plot <- function(model, title) {
  plot(model, 
       xlab = "", ylab = "", main = title, theme.bw = TRUE,
       Xdistr = "density", 
       cex.axis = 0.8, xlim = c(1, 5))
}

vote_engagement_plot <- make_interaction_plot(vote_engagement, "Δ (Vote for Han)")
thermo_engagement_plot <- make_interaction_plot(thermo_engagement, "Δ (Candidate Evaluations)")
china_engagement_plot <- make_interaction_plot(china_engagement, "Δ (Pro-PRC Index)")
combine_plot_interaction <- plot_grid(vote_engagement_plot, 
                                      thermo_engagement_plot, 
                                      china_engagement_plot, ncol = 3) # combine into one fig
combine_plot_interaction
ggsave(file.path(figures_dir, "Figure 5.pdf"), combine_plot_interaction, width = 12, height = 4, dpi = 300, device = my_pdf_device)


## Figure 6: Treatment Effects by Political Predispositions
# Function to process models for each outcome by partisan subgroups
make_subgroup_plot <- function(blue_itt, blue_att, 
                               indep_itt, indep_att, 
                               green_itt, green_att, title) {
  bind_rows(
    tidy(blue_itt) %>% filter(term != "(Intercept)") %>% mutate(group = "Pan-Blue", model = "ITT"),
    tidy(blue_att) %>% filter(term != "(Intercept)") %>% mutate(group = "Pan-Blue", model = "TOT"),
    tidy(indep_itt) %>% filter(term != "(Intercept)") %>% mutate(group = "Nonpartisan", model = "ITT"),
    tidy(indep_att) %>% filter(term != "(Intercept)") %>% mutate(group = "Nonpartisan", model = "TOT"),
    tidy(green_itt) %>% filter(term != "(Intercept)") %>% mutate(group = "Pan-Green", model = "ITT"),
    tidy(green_att) %>% filter(term != "(Intercept)") %>% mutate(group = "Pan-Green", model = "TOT")
  ) %>% mutate(title = title)
}

# Process all outcomes
subgroup.plot <- bind_rows(
  make_subgroup_plot(bluehan.itt, bluehan.att, indephan.itt, indephan.att, greenhan.itt, greenhan.att, "Δ Vote for Han"),
  make_subgroup_plot(bluefavhan.itt, bluefavhan.att, indepfavhan.itt, indepfavhan.att, greenfavhan.itt, greenfavhan.att, "Δ Candidate Evaluations"),
  make_subgroup_plot(bluechina.itt, bluechina.att, indepchina.itt, indepchina.att, greenchina.itt, greenchina.att, "Δ Pro-PRC Index")
)

# Set order
subgroup.plot$group <- factor(subgroup.plot$group, levels = c("Pan-Green", "Nonpartisan", "Pan-Blue"))
subgroup.plot$model <- factor(subgroup.plot$model, levels = c("TOT", "ITT"), 
                              labels = c("Treatment-on-the-Treated", "Intent-to-Treat"))
subgroup.plot$title <- factor(subgroup.plot$title, 
                              levels = c("Δ Vote for Han", "Δ Candidate Evaluations", "Δ Pro-PRC Index"))

ggplot(subgroup.plot, aes(x = group, y = estimate, shape = model)) +
  geom_hline(yintercept = 0, color = "gray60", linetype = "dashed") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.6),
                  size = 0.8) + 
  facet_grid(~ title, scales = "free_x") +
  scale_x_discrete(labels = c("Pan-Greens\n(PRC-skeptical)", "Nonpartisans", "Pan-Blues\n(PRC-friendly)")) +
  scale_shape_manual(values = c(1, 16)) +  
  theme_bw() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        legend.key.width = unit(2, "lines"),
        strip.background = element_blank(),
        strip.text = element_text(size = 14, face = "bold", margin = margin(b = 5)),
        axis.text = element_text(size = 12, color = "black"),
        axis.title.x = element_text(size = 12),
        panel.grid = element_blank()) +
  xlab("") + ylab("") + 
  coord_flip() +
  guides(shape = guide_legend(reverse = TRUE))

ggsave(file.path(figures_dir, "Figure 6.pdf"), width = 1100, height = 380, unit = "px", dpi = 96, device = my_pdf_device)



########################################################
############ Online Supplementary Materials ############ 
########################################################



## Figure SI-2: Sentiment Analysis of Treatment News ##

# -----------------------------------------------------------------------------
# Install Non-CRAN Packages (jiebaR)
# -----------------------------------------------------------------------------
# jiebaR (Chinese text segmentation) has been removed from CRAN.
# If you did not have this package, we must install it from GitHub using the 'remotes' package.


# Install 'jiebaRD' (Dictionary Data) FIRST
if (!requireNamespace("jiebaRD", quietly = TRUE)) {
  message("Installing 'jiebaRD' (dictionary data) from GitHub...")
  remotes::install_github("qinwf/jiebaRD", upgrade = "never", force = TRUE)
}

if (!requireNamespace("jiebaR", quietly = TRUE)) {
  message("Installing 'jiebaR' from GitHub...")
  remotes::install_github("qinwf/jiebaR", upgrade = "never", force = TRUE)
}

suppressPackageStartupMessages(library(jiebaR))


# --- I. Load Data ---
# Uses file.path(data_dir, ...) to ensure it works on any computer
ctnews <- read_csv(file.path(data_dir, "ctnews.csv"), show_col_types = FALSE)

# Load dictionaries
pos <- readLines(file.path(data_dir, "sentiment_positive.txt"), encoding = 'UTF-8', warn = FALSE)
neg <- readLines(file.path(data_dir, "sentiment_negative.txt"), encoding = 'UTF-8', warn = FALSE)
stopwords <- readLines(file.path(data_dir, "sentiment_stopword.txt"), encoding = 'UTF-8', warn = FALSE)

# --- II. Initialize Jieba & Custom Terms ---
# Note: jiebaR package must be installed
wk <- worker()

custom_terms <- c(
  '韓國瑜', '蔡英文', '反滲透法', '知識藍', '經濟藍', 
  '中聯辦', '國台辦', '年報', '波特王', '南台灣', 
  '高市府', '已經', '傷透了心')

# Add words one by one
walk(custom_terms, ~new_user_word(wk, .x))

# --- III. Processing Pipeline ---
ctnews_scored <- ctnews %>%
  mutate(
    date = mdy(date),
    # Calculate length
    length = nchar(content),
    # Tokenize
    words = map(content, ~segment(.x, wk)),
    # Remove Stopwords
    words = map(words, ~ .x[!.x %in% stopwords]),
    # Calculate Scores
    pos_weight = map_int(words, ~ sum(.x %in% pos)),
    neg_weight = map_int(words, ~ sum(.x %in% neg)),
    # Calculate Net Sentiment
    total = pos_weight - neg_weight,
    emotion = if_else(total >= 0, 'Pos', 'Neg'))

# --- IV. Visualization ---

# Count of Articles by Sentiment
p1 <- ggplot(ctnews_scored, aes(x = emotion)) +
  geom_bar(stat = "count", fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), vjust = -1) +
  labs(title = "Total Articles by Sentiment", x = "Sentiment", y = "Count") +
  theme_minimal()

# Sentiment by Publication Date
p2 <- ggplot(ctnews_scored, aes(x = date, fill = emotion)) +
  geom_bar(position = "fill") +
  labs(title = "Sentiment Proportion by Date", x = "Date", y = "Proportion") +
  scale_x_date(date_labels = "%m/%d", date_breaks = "1 day") + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Sentiment by Article Length
p3 <- ggplot(ctnews_scored, aes(x = length, fill = emotion)) +
  geom_density(alpha = 0.5) +
  labs(title = 'Sentiment by Article Length', x = 'Length', y = 'Density') +
  theme_minimal()

# Combine plots using patchwork (p1 | (p2 / p3))
figure_si2 <- p1 | (p2 / p3)

# Save the figure
ggsave(file.path(figures_dir, "Figure SI-2.pdf"), figure_si2, width = 12, height = 8)


## Figure SI-3: Distribution of Average Browsing Time Per Visit ##

website.df <-experiment.df %>%
  mutate(avgminute = avg_browsing/60,
         treatment = factor(treatment)) %>% 
  filter(conditions == "Placebo" | conditions == "Treatment") 

lm(avgminute ~ treatment, data = website.df) %>% summary()

ggplot(subset(website.df, avgminute<15), aes(x = avgminute, y = after_stat(density))) +  # Updated syntax
  geom_histogram(aes(fill = treatment), alpha = 0.6, binwidth = 0.3) +
  facet_grid(treatment ~ ., labeller = as_labeller(c("0" = "Placebo Website", "1" = "Treatment Website"))) +
  scale_x_continuous(name = "Average Time Spent Per Visit (Minutes)", breaks = seq(0, 15, by = 0.5)) +
  scale_y_continuous(name = "Density") +
  scale_fill_manual(values = c("0" = "#1f77b4", "1" = "#ff7f0e")) +  # Named values for clarity
  theme_bw() +
  guides(fill = "none")  

ggsave(file.path(figures_dir, "Figure SI-3.pdf"), width = 12.3, height = 5.68)

## Figure SI-4: Panel View of Browsing Pattern across Individual and Time ##

panel_tracking <- read_csv(file.path(data_dir, "panelvisiting.csv"),
                           col_types = cols(
                             gid = col_double(),
                             date = col_date(format = "%m/%d/%Y"), 
                             visited = col_double()))

panelview(date ~ visited, data = panel_tracking, 
          index = c("gid","date"), xlab = "Publication Dates", ylab = "Participants",
          color = c("white","red2"),
          legend.labs = c("Not Visited","Visited"),
          main = "Browsing Behavior on the Pro-Beijing News Website") 

ggsave(file.path(figures_dir, "Figure SI-4.pdf"), width = 10.5, height = 5.05)


## Figure SI-5: Predicted Probability of Minimum Compliance ##
# Function to fit model and create plot
make_pred.compliance_plot <- function(data, outcome_var, title = "Predicted Probabilities of Compliers by Partisanship") {
  
  # Fit model
  model <- glm(reformulate("spectrum", outcome_var), 
               family = binomial(link = "logit"), 
               data = data)
  # Get predictions
  pred <- predict(model, type = "response", se.fit = TRUE)
  # Create summary data frame
  summary_table <- data.frame(
    x = data$spectrum,
    Predicted_Probabilities = pred$fit,
    Lower_CI = pred$fit - 1.96 * pred$se.fit,
    Upper_CI = pred$fit + 1.96 * pred$se.fit)
  # Create plot
  ggplot(summary_table, aes(x = x, y = Predicted_Probabilities)) +
    geom_point(color = "darkblue", size = 4, shape = 21, fill = "lightblue") +
    geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI), 
                  width = 0.2, color = "darkred", linewidth = 1.2) +
    geom_line(aes(group = 1), color = "darkblue", linetype = "solid", linewidth = 1.2) +
    labs(title = "",
         subtitle = "With 95% Confidence Intervals (red bars)",
         x = "Political Predispositions",
         y = "Predicted Probability") +
    theme_minimal(base_size = 15) +
    theme(plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5))
}

# Subset data for Panel b of SI-5 and SI-6
experiment_wave2 <- subset(experiment.df, wave2 == 1)

# Figure SI-5: Minimum Compliance
plot_min_all <- make_pred.compliance_plot(experiment.df, "min_complier")
plot_min_completed <- make_pred.compliance_plot(experiment_wave2, "min_complier")
combine_min_compliance <- plot_grid(plot_min_all, plot_min_completed, ncol = 2)

ggsave(file.path(figures_dir, "Figure SI-5.pdf"), combine_min_compliance, width = 12, height = 4, dpi = 300)


# Figure SI-6: Full Compliance
plot_full_all <- make_pred.compliance_plot(experiment.df, "full_complier")
plot_full_completed <- make_pred.compliance_plot(experiment_wave2, "full_complier")
combine_full_compliance <- plot_grid(plot_full_all, plot_full_completed, ncol = 2)

ggsave(file.path(figures_dir, "Figure SI-6.pdf"), combine_full_compliance, width = 12, height = 4, dpi = 300)


## Figure SI-7: PRC Index Baseline Scores by Political Predispositions ##

# Function to create bar plot with error bars for PRC opinions by political groups
make_china_plot <- function(data, yvar, title) {
  base_china <- data %>% 
    group_by(spectrum) %>%
    do(tidy(lm_robust(reformulate("1", yvar), data = .))) %>%
    mutate(yval = estimate)
  
  ggplot(base_china, aes(x = spectrum, y = yval)) + 
    geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .1, linewidth = 1) +
    theme_bw() +
    xlab("") + ylab("Mean Score") + ggtitle(title)
}

fig_a7_1 <- make_china_plot(analysis, "china_thermo_w1", "Favorability toward China")
fig_a7_2 <- make_china_plot(analysis, "ch_econ_inf_w1", "Impact of PRC's Growth")
fig_a7_3 <- make_china_plot(analysis, "ch_threat_w1", "China Not a Threat")
fig_a7_4 <- make_china_plot(analysis, "ch_usa_power_w1", "China Won't Overtake the U.S.")
fig_a7_5 <- make_china_plot(analysis, "strait_trade_w1", "Trade with China")
fig_a7_6 <- make_china_plot(analysis, "xidada_w1", "Xi's Global Leadership")
fig_a7_7 <- make_china_plot(analysis, "hk_protest_w1", "Don't Support HK Protests")
fig_a7_8 <- make_china_plot(analysis, "hk_radical_w1", "Don't Support HK's Radical Protesters")

combine_base_china <- plot_grid(fig_a7_1, fig_a7_2, fig_a7_3, fig_a7_4, 
                                fig_a7_5, fig_a7_6, fig_a7_7, fig_a7_8, ncol = 4)

ggsave(file.path(figures_dir, "Figure SI-7.pdf"), combine_base_china, width = 14, height = 6, dpi = 300)



## Figure SI-8: Outcome Baseline Scores by Experimental Conditions ##

# Function to create bar plot with error bars for outcome's baseline score by conditions
make_baseline_plot <- function(data, yvar, title) {
  balanced_baseline <- data %>% 
    filter(conditions!= "Existing") %>%
    group_by(conditions) %>%
    do(tidy(lm_robust(reformulate("1", yvar), data = .))) %>%
    mutate(yval = estimate)
  
  ggplot(balanced_baseline, aes(x = conditions, y = yval)) + 
    geom_point(position = position_dodge(width = .7)) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
    theme_bw() +
    geom_text(aes(label = round(estimate, 2), y = estimate), 
              vjust = -1.5, hjust = -.1, size = 3.5, color = "black") +
    xlab("") + ylab("Mean Score") + ggtitle(title)
}

fig_a8_1 <- make_baseline_plot(experiment.df, "turnout_w1", "Turnout Intent")
fig_a8_2 <- make_baseline_plot(experiment.df, "vote_han_w1", "Vote Intent for Han")
fig_a8_3 <- make_baseline_plot(experiment.df, "vote_tsai_w1", "Vote Intent for Tsai")
fig_a8_4 <- make_baseline_plot(experiment.df, "vote_song_w1", "Vote Intent for Soong")
fig_a8_5 <- make_baseline_plot(experiment.df, "relative_thermo_w1", "Candidate Evaluations")
fig_a8_6 <- make_baseline_plot(experiment.df, "prc_index_w1", "Pro-PRC Index")
fig_a8_7 <- make_baseline_plot(analysis, "china_thermo_w1", "Favorability toward China")
fig_a8_8 <- make_baseline_plot(analysis, "ch_econ_inf_w1", "Impact of PRC's Growth")
fig_a8_9 <- make_baseline_plot(analysis, "ch_threat_w1", "China Not a Threat")
fig_a8_10 <- make_baseline_plot(analysis, "ch_usa_power_w1", "China Won't Overtake the U.S.")
fig_a8_11 <- make_baseline_plot(analysis, "strait_trade_w1", "Trade with China")
fig_a8_12 <- make_baseline_plot(analysis, "xidada_w1", "Xi's Global Leadership")
fig_a8_13 <- make_baseline_plot(analysis, "hk_protest_w1", "Don't Support HK Protests")
fig_a8_14 <- make_baseline_plot(analysis, "hk_radical_w1", "Don't Support HK's Radical Protesters")

combine_balanced_base <- 
  plot_grid(fig_a8_1, fig_a8_2, fig_a8_3, fig_a8_4, fig_a8_5, 
            fig_a8_6, fig_a8_7, fig_a8_8, fig_a8_9, fig_a8_10,
            fig_a8_11, fig_a8_12, fig_a8_13, fig_a8_14, ncol = 3)
ggsave(file.path(figures_dir, "Figure SI-8.pdf"), combine_balanced_base, width = 18, height = 10, dpi = 300)


## Figure SI-9: Predicted Endline Survey Participation ##

# Break age into groups for conservative estimates
experiment.df <- experiment.df %>%
  mutate(age_group = case_when(
    age >= 20 & age <= 24 ~ 1,
    age >= 25 & age <= 34 ~ 2,
    age >= 35 & age <= 44 ~ 3,
    age >= 45 & age <= 54 ~ 4,
    age >= 55 ~ 5),
    age_group = factor(age_group, 
                       levels = 1:5, 
                       labels = c("20-24", "25-34", "35-44", "45-54", "55+")))

attrition_interaction <- glm(wave2 ~ 
                               conditions * age_group + 
                               conditions * female +
                               conditions * special_municipality +
                               conditions * edu + 
                               conditions * full_time +
                               conditions * married + 
                               conditions * income_cat +
                               conditions * spectrum + 
                               conditions * president_2016_voted +
                               conditions * president_2012_voted + 
                               conditions * father_islander +
                               conditions * taoism + 
                               conditions * nwpaper_pref_apple +
                               conditions * sanlih_show,
                             family = binomial(link = "probit"),
                             data=experiment.df)

attrition_interaction.tidy <- tidy(attrition_interaction)

dwplot(attrition_interaction.tidy,
       dot_args = list(color="black"),
       whisker_args = list(color="black", size = 0.6)) %>% 
  relabel_predictors(c(conditionsPlacebo = "Placebo",
                       conditionsTreatment = "Treatment",
                       "age_group25-34"= "Age 25-34",
                       "conditionsPlacebo:age_group25-34" = "Placebo x Age 25-34",
                       "conditionsTreatment:age_group25-34" = "Treatment x Age 25-34",
                       "age_group35-44" = "Age 35-44",
                       "conditionsPlacebo:age_group35-44" = "Placebo x Age 35-44",
                       "conditionsTreatment:age_group35-44" = "Treatment x Age 35-44",
                       "age_group45-54" = "Age 45-54",
                       "conditionsPlacebo:age_group45-54" = "Placebo x Age 45-54",
                       "conditionsTreatment:age_group45-54" = "Treatment x Age 45-54",
                       "age_group55+" = "Age 55+",
                       "conditionsPlacebo:age_group55+" = "Placebo x Age 55+",
                       "conditionsTreatment:age_group55+" = "Treatment x Age 55+",
                       female = "Female",
                       "conditionsPlacebo:female" = "Placebo x Female",
                       "conditionsTreatment:female" = "Treatment x Female",
                       special_municipality = "Residence",
                       "conditionsPlacebo:special_municipality" = "Placebo x Residence",
                       "conditionsTreatment:special_municipality" = "Treatment x Residence",
                       edu = "Education",
                       "conditionsPlacebo:edu" = "Placebo x Education",
                       "conditionsTreatment:edu" = "Treatment x Education",
                       full_time = "Full-time Job",
                       "conditionsPlacebo:full_time" = "Placebo x Full-time Job",
                       "conditionsTreatment:full_time" = "Treatment x Full-time Job",
                       married = "Married",
                       "conditionsPlacebo:married" = "Placebo x Married",
                       "conditionsTreatment:married" = "Treatment x Married",
                       income_cat = "Income",
                       "conditionsPlacebo:income_cat" = "Placebo x Income",
                       "conditionsTreatment:income_cat" = "Treatment x Income",
                       spectrumNonpartisan = "Nonpartisan Voter",
                       "conditionsPlacebo:spectrumNonpartisan" = "Placebo x Nonpartisan Voter",
                       "conditionsTreatment:spectrumNonpartisan" = "Treatment x Nonpartisan Voter",
                       "spectrumPan-Green" = "Pan-green Voter",
                       "conditionsPlacebo:spectrumPan-Green" = "Placebo x Pan-green Voter",
                       "conditionsTreatment:spectrumPan-Green" = "Treatment x Pan-green Voter",
                       president_2016_voted = "Vote in 2016",
                       "conditionsPlacebo:president_2016_voted" = "Placebo x Vote in 2016",
                       "conditionsTreatment:president_2016_voted" = "Treatment x Vote in 2016",
                       president_2012_voted = "Vote in 2012",
                       "conditionsPlacebo:president_2012_voted" = "Placebo x Vote in 2012",
                       "conditionsTreatment:president_2012_voted" = "Treatment x Vote in 2012",
                       father_islander = "Minnan",
                       "conditionsPlacebo:father_islander" = "Placebo x Minnan",
                       "conditionsTreatment:father_islander" = "Treatment x Minnan",
                       taoism = "Taoism",
                       "conditionsPlacebo:taoism" = "Placebo x Taoism",
                       "conditionsTreatment:taoism" = "Treatment x Taoism",
                       nwpaper_pref_apple = "Newspaper Preference",
                       "conditionsPlacebo:nwpaper_pref_apple" = "Placebo x Newspaper Preference",
                       "conditionsTreatment:nwpaper_pref_apple" = "Treatment x Newspaper Preference",
                       sanlih_show = "TV Program Preference",
                       "conditionsPlacebo:sanlih_show" = "Placebo x TV Program Preference",
                       "conditionsTreatment:sanlih_show" = "Treatment x TV Program Preference")) +
  xlab("Coefficients") + 
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme(legend.position = "none",
        axis.text=element_text(size=10),
        axis.title = element_text(face="bold")) + theme_bw()
ggsave(file.path(figures_dir, "Figure SI-9.pdf"), width = 8.5, height = 9.65)


## Figure SI-10: Placebo-Control Comparison ##

control.placebo.df <- experiment.df %>%
  filter(conditions == 'Placebo' | conditions == 'Control')

# Define outcomes and labels
outcomes <- c("turnout_change", "vote_han_change", "vote_tsai_change", "vote_song_change", 
              "relative_thermo_change", "prc_index_change",
              "china_thermo_change", "ch_econ_inf_change", "ch_threat_change",
              "ch_usa_power_change", "strait_trade_change", "xidada_change",
              "hk_protest_change", "hk_radical_change")

labels <- c("Δ (Turnout)", "Δ (Vote for Han)", "Δ (Vote for Tsai)", "Δ (Vote for Soong)",
            "Δ (Candidate Evaluation)", "Δ (Pro-China Index)",
            "Δ (Favor. China)", "Δ (Chinese Economy)", "Δ (China Threat)",
            "Δ (Overtake U.S.)", "Δ (Trade w/ China)", "Δ (Xi Global Leadership)",
            "Δ (HK Protest)", "Δ (Radical Behavior)")

# Run models and combine results
results <- map_dfr(outcomes, function(y) {
  lm_robust(reformulate("placebo", y), data = control.placebo.df) %>%
    tidy() %>%
    filter(term != "(Intercept)") %>%
    mutate(outcome = y)
})

# Set factor levels
results$outcome <- factor(results$outcome, levels = outcomes)

# Plot
ggplot(results, aes(x = outcome, y = estimate)) + 
  geom_point(size = 4, color = "#d62728") + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, color = "#d62728") + 
  geom_hline(yintercept = 0, color = "gray", linetype = 'dashed') +
  theme_bw() +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 10, colour = "black"),
        axis.text.y = element_text(size = 12, colour = "black"),
        axis.title = element_text(size = 14),
        panel.grid = element_blank(),
        plot.margin = unit(c(1, 1, 1, 1), "cm")) +
  xlab("Outcomes") + ylab("Estimates") +
  scale_x_discrete(labels = labels)

ggsave(file.path(figures_dir, "Figure SI-10.pdf"), width = 22, height = 8, device = my_pdf_device)


## Figure SI-11: Pre-election Polling of the General Election ##

# Read polling data
polling <- read_csv(file.path(data_dir, "polling.csv"), # Updated Input Path
                    show_col_types = FALSE) %>%
  mutate(
    across(c(han, tsai, song), ~ as.numeric(gsub("%", "", .)) / 100),
    date = as.Date(date, format = "%m/%d/%Y"))

polling.long <- polling %>%
  pivot_longer(cols = c(han, tsai, song), 
               names_to = "candidate", 
               values_to = "results") %>%
  mutate(candidate = factor(candidate, 
                            levels = c("han", "tsai", "song"),
                            labels = c("Han", "Tsai", "Soong"))) %>%
  filter(date != as.Date("2020-01-11"))

# Date breaks and labels
date_breaks <- as.Date(c('2019-10-18', '2019-10-28', '2019-11-08', '2019-11-18', 
                         '2019-11-28', '2019-12-08', '2019-12-18', '2019-12-28', '2020-01-06'))
date_labels <- c('2019-10-18', '2019-10-28', '2019-11-08', '2019-11-18', 
                 '2019-11-28', '2019-12-08', '2019-12-18', '2019-12-28', 'Election Day\n(2020-01-11)')

# Annotations
election_results <- data.frame(
  x = as.Date("2020-01-06"),
  y = c(.5613, .3761, .0526),
  color = c("seagreen3", "steelblue1", "tan1")
)

vote_labels <- data.frame(
  x = as.Date("2020-01-12"),
  y = c(.56, .375, .06),
  label = c("Tsai's\nVote Share", "Han's\nVote Share", "Soong's\nVote Share")
)

# Plot
ggplot(polling.long, aes(x = date, y = results, color = candidate)) +
  geom_point(aes(shape = candidate), size = 2) +
  geom_smooth(span = .4, se = TRUE) +
  geom_vline(xintercept = as.numeric(as.Date("2020-01-06")), linetype = 2) +
  geom_point(data = election_results, aes(x = x, y = y), 
             color = election_results$color, size = 4, shape = 10, inherit.aes = FALSE) +
  geom_label(data = vote_labels, aes(x = x, y = y, label = label), inherit.aes = FALSE) +
  scale_color_manual(name = "Candidates", values = c("steelblue1", "seagreen3", "tan1")) +
  scale_x_date(breaks = date_breaks, labels = date_labels,
               limits = as.Date(c('2019-10-15', '2020-01-13'))) +
  scale_y_continuous(breaks = scales::breaks_extended(8), labels = scales::label_percent()) +
  guides(shape = "none") +
  labs(x = "Polling Dates", y = "Percentage", title = "") +
  theme_bw()

ggsave(file.path(figures_dir, "Figure SI-11.pdf"), width = 12, height = 6)


## Figure SI-12: Treatment Effects on Separate PRC Index Items ##


# Function to process models for each outcome
make_china_item_plot <- function(existing, itt, min, title) {
  bind_rows(
    tidy(existing) %>% mutate(model = "Existing"),
    tidy(itt) %>% filter(term != "(Intercept)") %>% mutate(model = "ITT"),
    tidy(min) %>% filter(term != "(Intercept)") %>% mutate(model = "Minimal")
  ) %>% mutate(title = title)
}

# Fit models for each outcome
# A: Favorability of China
favorchina.existing <- lm_robust(china_thermo_change ~ 1, data = experiment.df)
favorchina.itt <- lm_robust(china_thermo_change ~ treatment, data = experiment.df)
favorchina.min <- iv_robust(china_thermo_change ~ min_complier | treatment, data = experiment.df)

# B: China Threat
threat.existing <- lm_robust(ch_threat_change ~ 1, data = experiment.df)
threat.itt <- lm_robust(ch_threat_change ~ treatment, data = experiment.df)
threat.min <- iv_robust(ch_threat_change ~ min_complier | treatment, data = experiment.df)

# C: Trade with China
trade.existing <- lm_robust(strait_trade_change ~ 1, data = experiment.df)
trade.itt <- lm_robust(strait_trade_change ~ treatment, data = experiment.df)
trade.min <- iv_robust(strait_trade_change ~ min_complier | treatment, data = experiment.df)

# D: Support for HK
hk.existing <- lm_robust(hk_protest_change ~ 1, data = experiment.df)
hk.itt <- lm_robust(hk_protest_change ~ treatment, data = experiment.df)
hk.min <- iv_robust(hk_protest_change ~ min_complier | treatment, data = experiment.df)

# E: Radical Behavior of HK Protests
radical.existing <- lm_robust(hk_radical_change ~ 1, data = experiment.df)
radical.itt <- lm_robust(hk_radical_change ~ treatment, data = experiment.df)
radical.min <- iv_robust(hk_radical_change ~ min_complier | treatment, data = experiment.df)

# F: Influence of Chinese economy
inf.existing <- lm_robust(ch_econ_inf_change ~ 1, data = experiment.df)
inf.itt <- lm_robust(ch_econ_inf_change ~ treatment, data = experiment.df)
inf.min <- iv_robust(ch_econ_inf_change ~ min_complier | treatment, data = experiment.df)

# G: China overtakes US
overtake.existing <- lm_robust(ch_usa_power_change ~ 1, data = experiment.df)
overtake.itt <- lm_robust(ch_usa_power_change ~ treatment, data = experiment.df)
overtake.min <- iv_robust(ch_usa_power_change ~ min_complier | treatment, data = experiment.df)

# H: Xi's global leadership
xidada.existing <- lm_robust(xidada_change ~ 1, data = experiment.df)
xidada.itt <- lm_robust(xidada_change ~ treatment, data = experiment.df)
xidada.min <- iv_robust(xidada_change ~ min_complier | treatment, data = experiment.df)

# Combine all results
china.sep <- bind_rows(
  make_china_item_plot(favorchina.existing, favorchina.itt, favorchina.min, "Favorability toward China"),
  make_china_item_plot(inf.existing, inf.itt, inf.min, "Influence of Chinese Economy"),
  make_china_item_plot(threat.existing, threat.itt, threat.min, "China Threat"),
  make_china_item_plot(overtake.existing, overtake.itt, overtake.min, "China Overtakes the U.S."),
  make_china_item_plot(trade.existing, trade.itt, trade.min, "Trade with China"),
  make_china_item_plot(xidada.existing, xidada.itt, xidada.min, "Xi's Global Leadership"),
  make_china_item_plot(hk.existing, hk.itt, hk.min, "Support for HK Protests"),
  make_china_item_plot(radical.existing, radical.itt, radical.min, "Radical Behavior in HK Protests")
)

# Set factor levels
china.sep$title <- factor(china.sep$title, 
                          levels = c("Favorability toward China", 
                                     "Influence of Chinese Economy",
                                     "China Threat",
                                     "China Overtakes the U.S.",
                                     "Trade with China",
                                     "Xi's Global Leadership",
                                     "Support for HK Protests",
                                     "Radical Behavior in HK Protests"))

china.sep$model <- factor(china.sep$model, levels = c("Existing", "Minimal", "ITT"))

# Plot
ggplot(china.sep, aes(x = model, y = estimate, fill = model)) + 
  geom_point(size = 3.5, shape = 21, color = "black") + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, linewidth = 0.8) + 
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  scale_x_discrete(labels = c("Existing", "TOT", "ITT")) +
  facet_wrap(~ title, scales = "free", nrow = 2, ncol = 4) +
  scale_fill_manual(values = c("#1f77b4", "#2ca02c", "#d62728")) +
  theme_bw() +
  theme(legend.position = "none",
        strip.background = element_rect(fill = "white"),
        strip.text = element_text(size = 14, face = "bold", margin = margin(b = 5)),
        axis.text.x = element_text(size = 10, colour = "black", angle = 45, hjust = 1),
        axis.text.y = element_text(size = 12, colour = "black"),
        panel.spacing = unit(1, "lines"),
        panel.grid = element_blank(),
        plot.margin = margin(10, 10, 10, 10)) +
  xlab("") + ylab("Estimate") + 
  coord_flip()

ggsave(file.path(figures_dir, "Figure SI-12.pdf"), width = 16, height = 6)



## Figure SI-13: Cognitive and Affective Reactions to Treatment Site News ##

# Subset data ((PRC-skeptical vs. PRC-friendly in treatment site)
experiment.treatment <- subset(experiment.df, conditions == "Treatment" & spectrum != "Nonpartisan")

# Function to create bar plot by partisanship 
make_mechanism_plot <- function(data, yvar, title) {
  
  summary_df <- data %>% 
    group_by(spectrum) %>%
    do(tidy(lm_robust(reformulate("1", yvar), data = .))) %>%
    mutate(yval = estimate)
  
  ggplot(summary_df, aes(x = spectrum, y = yval)) + 
    geom_bar(stat = "identity", fill = c('#2874A6', '#239B56')) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
    geom_text(aes(label = round(yval, 2), y = conf.high), vjust = 5.5, size = 3, color = "white") +
    xlab("") + ylab("Mean Score") + ggtitle(title) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5),
          panel.grid = element_blank())
}


# Create plots
plot_quality <- make_mechanism_plot(experiment.treatment, "website_quality", "Quality")
plot_strength <- make_mechanism_plot(experiment.treatment, "website_strength", "Argument Strength")
plot_negative <- make_mechanism_plot(experiment.treatment, "negative_emotion", "Negative Feeling")
plot_positive <- make_mechanism_plot(experiment.treatment, "positive_emotion", "Positive Feeling")
plot_reactance <- make_mechanism_plot(experiment.treatment, "reactance", "Reactance")

# Combined plots
website_mechanism <-(plot_quality + plot_strength) / (plot_positive + plot_negative + plot_reactance) +
  plot_annotation(title = "Treatment Website")
website_mechanism
ggsave(file.path(figures_dir, "Figure SI-13.pdf"), website_mechanism, width = 8.5, height = 6.5)


## Figure SI-14: Cognitive and Affective Reactions to Placebo Site News ##

# Subset data ((PRC-skeptical vs. PRC-friendly in placebo site)
experiment.placebo <- subset(experiment.df, conditions == "Placebo" & spectrum != "Nonpartisan")

plot_quality_placebo <- make_mechanism_plot(experiment.placebo, "website_quality", "Quality")
plot_strength_placebo <- make_mechanism_plot(experiment.placebo, "website_strength", "Argument Strength")
plot_negative_placebo <- make_mechanism_plot(experiment.placebo, "negative_emotion", "Negative Feeling")
plot_positive_placebo <- make_mechanism_plot(experiment.placebo, "positive_emotion", "Positive Feeling")
plot_reactance_placebo <- make_mechanism_plot(experiment.placebo, "reactance", "Reactance")

website_mechanism_placebo <-(plot_quality_placebo + plot_strength_placebo) / 
  (plot_positive_placebo + plot_negative_placebo + plot_reactance_placebo) +
  plot_annotation(title = "Placebo Website")
website_mechanism_placebo
ggsave(file.path(figures_dir, "Figure SI-14.pdf"), website_mechanism_placebo, width = 8.5, height = 6.5)


## Figure SI-15: Treatment-Placebo Comparison ##

tp_vote.itt <- lm_robust(vote_han_change ~ treatment, data = subset(experiment.df, conditions!="Control"))
tp_favhan.itt <- lm_robust(relative_thermo_change ~ treatment, data = subset(experiment.df, conditions!="Control"))
tp_china.itt <- lm_robust(prc_index_change ~ treatment, data = subset(experiment.df, conditions!="Control"))

# Combine results
tp_full.plot <- bind_rows(
  tidy(tp_vote.itt) %>% filter(term != "(Intercept)") %>% mutate(outcome = "Δ Vote for Han"),
  tidy(tp_favhan.itt) %>% filter(term != "(Intercept)") %>% mutate(outcome = "Δ Candidate Evaluations"),
  tidy(tp_china.itt) %>% filter(term != "(Intercept)") %>% mutate(outcome = "Δ Pro-PRC Index")
)

# Set factor levels
tp_full.plot$outcome <- factor(tp_full.plot$outcome, 
                               levels = c("Δ Pro-PRC Index",
                                          "Δ Candidate Evaluations",
                                          "Δ Vote for Han"))

ggplot(tp_full.plot, aes(x = estimate, y = outcome)) + 
  geom_point(size = 3) + 
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0) + 
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray40") +
  theme_bw() + 
  geom_text(aes(label = round(estimate, 3)), 
            nudge_y = 0.25, 
            size = 3.5, 
            hjust = 0) +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        plot.margin = margin(10, 10, 10, 10)) +
  xlab("Treatment Effects (Treatment - Placebo)") + ylab("") 

ggsave(file.path(figures_dir, "Figure SI-15.pdf"), width = 8, height = 4.5, device = my_pdf_device)


## Figure SI-16: List Experiment ##
list.exp <- experiment.df %>%
  filter(wave2 == 1) %>%
  mutate(list_exp_group = case_when(
    list_exp_group == 0 ~ 'Control\n(N=483)',
    list_exp_group == 1 ~ 'Treatment\n(N=466)'))

summary <- list.exp %>% 
  group_by(list_exp_group) %>%
  do(tidy(lm_robust(itemcount ~ 1, data = .))) %>%
  mutate(itemcount = estimate)

ggplot(list.exp, aes(list_exp_group, itemcount)) +
  geom_point(position = position_jitter(width = .25, height = .25, seed = 0316), 
             alpha = 0.1, stroke = 0) +
  geom_pointrange(data = summary, 
                  aes(y = itemcount, ymin = conf.low, ymax = conf.high),
                  size = 0.5) +
  theme_bw() +
  theme(strip.background = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "none",
        panel.grid = element_blank()) +
  xlab("Item Counts") + ylim(1, 5)

ggsave(file.path(figures_dir, "Figure SI-16.pdf"), width= 4 , height=6, dpi = 300)


## Figure SI-17: Bounded Estimates of Treatment Effects ##

#Impute no change to attriters
bounded.df <- experiment.df %>%
  mutate(across(c(vote_han_change, relative_thermo_change, prc_index_change), 
                ~ replace_na(., 0)))

make_bounded_plot <- function(outcome, data, title) {
  itt <- lm_robust(reformulate("treatment", outcome), data = data)
  tot <- iv_robust(as.formula(paste(outcome, "~ full_complier | treatment")), data = data)
  
  bind_rows(
    tidy(itt) %>% filter(term != "(Intercept)") %>% mutate(model = "ITT"),
    tidy(tot) %>% filter(term != "(Intercept)") %>% mutate(model = "TOT")
  ) %>% mutate(title = title)
}

# For each outcome
bounded.plot <- bind_rows(
  make_bounded_plot("vote_han_change", bounded.df, "Δ (Vote for Han)"),
  make_bounded_plot("relative_thermo_change", bounded.df, "Δ (Candidate Evaluations)"),
  make_bounded_plot("prc_index_change", bounded.df, "Δ (Pro-PRC Index)")
)

# Set order
bounded.plot$title <- factor(bounded.plot$title, 
                             levels = c("Δ (Vote for Han)", 
                                        "Δ (Candidate Evaluations)",
                                        "Δ (Pro-PRC Index)"))

bounded.plot$model <- factor(bounded.plot$model, levels = c("TOT", "ITT"))

# Plot
ggplot(bounded.plot, aes(x = model, y = estimate, fill = model)) + 
  geom_point(size = 4, shape = 21, color = "black") + 
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, linewidth = 1) + 
  geom_hline(yintercept = 0, linetype = 2, color = "gray40") +
  scale_x_discrete(labels = c("Treatment-on-the-Treated", "Intent-to-Treat")) +
  facet_grid(. ~ title, scales = "free") +
  scale_fill_manual(values = c("#1f77b4", "#d62728")) +
  theme_bw() +
  theme(legend.position = "none",
        strip.background = element_blank()) +
  xlab(NULL) + ylab(NULL) +
  coord_flip()
ggsave(file.path(figures_dir, "Figure SI-17.pdf"), width = 12, height = 5, device = my_pdf_device)


## Figure SI-18: Treatment Effects w/ Inverse Probability Weighting ##

# Get predicted probabilities and compute IPW weights
experiment.df <- experiment.df %>%
  mutate(p_stayed = predict(attrition_interaction, type = "response"),
         ipw = 1 / p_stayed)

# Create survey design
design <- svydesign(ids = ~1, weights = ~ipw, data = experiment.df)

# Re-define outcomes and labels
outcomes <- c("vote_han_change", "relative_thermo_change", "prc_index_change")
labels <- c("Δ Vote for Han", "Δ Candidate Evaluations", "Δ Pro-PRC Index")

# Estimate IPW models and combine results
ipw_results <- bind_rows(
  tidy(svyglm(vote_han_change ~ treatment, design = design), conf.int = TRUE) %>%
    filter(grepl("treatment", term)) %>% mutate(outcome = "Δ Vote for Han"),
  tidy(svyglm(relative_thermo_change ~ treatment, design = design), conf.int = TRUE) %>%
    filter(grepl("treatment", term)) %>% mutate(outcome = "Δ Candidate Evaluations"),
  tidy(svyglm(prc_index_change ~ treatment, design = design), conf.int = TRUE) %>%
    filter(grepl("treatment", term)) %>% mutate(outcome = "Δ Pro-PRC Index")
)

# Set factor levels
ipw_results$outcome <- factor(ipw_results$outcome, 
                              levels = c("Δ Pro-PRC Index",
                                         "Δ Candidate Evaluations",
                                         "Δ Vote for Han"))

# Plot
ggplot(ipw_results, aes(x = estimate, y = outcome)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray40") +
  geom_text(aes(label = round(estimate, 3)), nudge_y = 0.25, size = 3.5, hjust = 0) +
  theme_bw() +
  theme(axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        panel.grid = element_blank()) +
  labs(x = "IPW-adjusted Treatment Effects", y = "")

# 635 x 355

# diagnosis of IPW distribution
summary(experiment.df$ipw)
ggsave(file.path(figures_dir, "Figure SI-18.pdf"), width = 8, height = 4.5, device = my_pdf_device)



## Figure SI-19: Distribution of IPW Weights ##
pdf(file.path(figures_dir, "Figure SI-19.pdf"), width = 7, height = 4)
hist(experiment.df$ipw, breaks = 30, main = "Distribution of IPW Weights", xlab = "Weight")
dev.off()



## Figure SI-20: Placebo Candidate ##
vote.soong.existing <- lm_robust(vote_song_change ~ 1, data = existing.df)
vote.soong.itt <- lm_robust(vote_song_change ~ treatment, data = experiment.df)
vote.soong.min <- iv_robust(vote_song_change ~ min_complier | treatment, data = experiment.df)
vote.soong.full <- iv_robust(vote_song_change ~ full_complier | treatment, data = experiment.df)

placebo_candidate <- make_main_plot(vote.soong.existing, vote.soong.itt, vote.soong.min, vote.soong.full, "Δ (Vote for Soong)")

ggplot(placebo_candidate, aes(x = model, y = estimate)) + 
  geom_hline(yintercept = 0, color = "gray40", linetype = "dashed") +
  geom_pointrange(aes(ymin = conf.low,
                      ymax = conf.high)) + 
  scale_x_discrete(labels = c("Existing CT\nConsumers",
                              "Treatment-on-the\nTreated (Full)",
                              "Treatment-on-the\nTreated (Minimal)",
                              "Intent-to-Treat")) +
  facet_grid(. ~ title, scales = "free") +
  coord_flip() + theme_bw() + 
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 14, face = "bold", margin = margin(b = 5)),
        axis.text = element_text(size = 12, color = "black"),
        axis.title = element_blank(),
        panel.grid = element_blank())

ggsave(file.path(figures_dir, "Figure SI-20.pdf"), width = 7.3, height = 4.35, device = my_pdf_device) #7.27 x 4.35



##Figure SI-21: Exploratory Analysis (Turnout) ##

turnout_all <- lm_lin(turnout_change ~ treatment, 
                      covariates = ~ pan_green + pan_blue + age + female + special_municipality +
                        edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                        father_islander + taoism + nwpaper_pref_apple + sanlih_show, 
                      data = experiment.df) 

# Effects on China-friendly subjects
turnout_blue <- lm_lin(turnout_change ~ treatment, 
                       covariates = ~ age + female + special_municipality +
                         edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                         father_islander + taoism + nwpaper_pref_apple + sanlih_show, 
                       data = subset(experiment.df, spectrum=="Pan-Blue"))
# Effects on nonpartisan subjects
turnout_nonpartisan <- lm_lin(turnout_change ~ treatment, 
                              covariates = ~ age + female + special_municipality +
                                edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                                father_islander + taoism + nwpaper_pref_apple + sanlih_show, 
                              data = subset(experiment.df, spectrum=="Nonpartisan")) 
# Effects on China-skeptical subjects
turnout_green <- lm_lin(turnout_change ~ treatment, 
                        covariates = ~ age + female + special_municipality +
                          edu + full_time + married + income_cat + president_2016_voted + president_2012_voted +
                          father_islander + taoism + nwpaper_pref_apple + sanlih_show, 
                        data = subset(experiment.df, spectrum=="Pan-Green")) 

turnout_effect <- bind_rows(
  "All Participants" = tidy(turnout_all),
  "PRC-Friendly"     = tidy(turnout_blue),
  "Nonpartisan"      = tidy(turnout_nonpartisan),
  "PRC-Skeptical"    = tidy(turnout_green),
  .id = "sample") %>%
  filter(term == "treatment") %>%
  select(sample, estimate, conf.low, conf.high) %>%
  mutate(sample = factor(sample, levels = c(
    "PRC-Skeptical", "Nonpartisan", "PRC-Friendly", "All Participants")))

# Plot
ggplot(turnout_effect, aes(x = estimate, y = sample)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
  labs(x = "Treatment Effects", y = NULL) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12))

ggsave(file.path(figures_dir, "Figure SI-21.pdf"), width = 635, height = 355, units = "px",, dpi = 96)



## Figure SI-22: Exploratory Analysis (Subgroup heterogeneity) ##

# One-time setup for subgroup identifiers
experiment.df <- experiment.df %>%
  mutate(tondu = case_when(
                 strait_relation >= 5 ~ 0,
                 strait_relation %in% c(2, 4) ~ 1,
                 strait_relation %in% c(1, 3) ~ 2),
        college = if_else(edu >= 5, 1, 0),
        taiwanese = if_else(identity ==1, 1, 0))

# Define the plot function
make_exp_subgroup_plot <- function(var_name, label_map, data = experiment.df) {
  
  # A. Run Model dynamically
  f <- as.formula(paste("vote_han_change ~ factor(treatment) * factor(", var_name, ")"))
  model <- lm_robust(f, data = data)
  
  # B. Calculate Margins & Clean Data
  margins <- ggpredict(model, terms = c("treatment", var_name)) %>%
    as.data.frame() %>%
    select(x, group, predicted, std.error) %>%
    pivot_wider(
      names_from = x, 
      values_from = c(predicted, std.error),
      names_sep = "_T"
    ) %>%
    mutate(
      treatment_effect = predicted_T1 - predicted_T0,
      se_diff = sqrt(std.error_T1^2 + std.error_T0^2), 
      ci_lower = treatment_effect - 1.96 * se_diff,
      ci_upper = treatment_effect + 1.96 * se_diff
    )
  
  # C. Plot
  p <- ggplot(margins, aes(x = group, y = treatment_effect)) +
    geom_hline(yintercept = 0, color = "gray40", linetype = "dashed") +
    geom_point(size = 3, color = "black") +
    geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0, color = "black", linewidth = 0.7) +
    scale_x_discrete(labels = label_map) +
    labs(x = NULL, y = NULL, title = NULL) +
    theme_bw() +
    theme(
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank()
    )
  return(p)
}

#Plot
# National Identity
make_exp_subgroup_plot(
  var_name = "taiwanese",
  label_map = c("0" = "Dual/Chinese", "1" = "Taiwanese"))

ggsave(file.path(figures_dir, "Figure SI-22A.pdf"), width = 658, height = 341, units = "px",, dpi = 96)

# Tondu
make_exp_subgroup_plot(
  var_name = "tondu",
  label_map = c("0" = "Status-Quo", "1" = "Independence", "2" = "Unification"))

ggsave(file.path(figures_dir, "Figure SI-22B.pdf"), width = 658, height = 341, units = "px",, dpi = 96)

# Ethnic Background
make_exp_subgroup_plot(
  var_name = "father_islander",
  label_map = c("0" = "People Outside the Province", "1" = "People of the Province"))

ggsave(file.path(figures_dir, "Figure SI-22C.pdf"), width = 658, height = 341, units = "px",, dpi = 96)


# Urban Residence
make_exp_subgroup_plot(
  var_name = "special_municipality",
  label_map = c("0" = "Living outside Special Municipalities", "1" = "Living in Special Municipalities"))

ggsave(file.path(figures_dir, "Figure SI-22D.pdf"), width = 658, height = 341, units = "px",, dpi = 96)

# College education
make_exp_subgroup_plot(
  var_name = "college",
  label_map = c("0" = "Without College Degree", "1" = "With College Degree"))

ggsave(file.path(figures_dir, "Figure SI-22E.pdf"), width = 658, height = 341, units = "px",, dpi = 96)



cat("Figures successfully reproduced. End of figure-replication script.\n")
